#################################
##' Code for:
##' Title: "Effective climate clubs require ambition, leverage, and insulation"
##' Author: Sam S. Rowan
##' Contact: sam.rowan [at] concordia.ca
##' Date: 2023-10-19
#################################
##' This script:
##' Baseline model
#################################


#### Initialize environment ####

## Packages
library(tidyverse)
library(truncnorm)
library(knitr)

## Working directory
# setwd("/Users/srowa/Dropbox/Projects/climate_clubs_trade")
# -- Set working directory to replication folder

## Define the logistic function
logistic_function <- function(x, k, a) {
  return(1 / (1 + exp(-k * (x - a))))
}

## Plot settings for later
theme_clubs <-  theme(legend.position = "bottom", 
                      panel.background = element_rect(fill=NA),
                      panel.grid.major.x = element_blank(),
                      panel.grid.major.y = element_blank(),
                      panel.grid.minor.y =  element_blank(),
                      panel.grid.minor.x =  element_blank(),
                      panel.ontop = FALSE,
                      axis.line = element_line(linewidth = 1, colour = "black"),
                      axis.text = element_text(size = 16),
                      legend.text = element_text(size = 12),
                      legend.title = element_text(size = 16),
                      axis.title = element_text(size = 16))

#### Define climate model parameters ####

## Define number of states
num_states <- 150


## Draw their climate ideal points
set.seed(202109)
climate_ideal <- cbind.data.frame(
  climate_ideal = runif(n = num_states, min = 0, max = 1)) %>% 
  arrange(climate_ideal) %>% 
  mutate(id = seq(from = 1, to = num_states, by = 1))
# climate_ideal


## Draw dyadic trade flows as a function of ideal points
climate_trade_data <- cbind.data.frame(
  # Add IDs
  id_a = rep(climate_ideal$id, times = num_states),
  id_b = rep(climate_ideal$id, each = num_states),
  # Bind climate ideal points into dyadic setup
  climate_ideal_a = rep(climate_ideal$climate_ideal, times = num_states),
  climate_ideal_b = rep(climate_ideal$climate_ideal, each = num_states)) %>%
  # Measure dyadic ideal point distance
  mutate(ideal_point_distance = abs(climate_ideal_a - climate_ideal_b)) %>%
  # Draw trade data based on ideal point distance
  mutate(trade_independent = rtruncnorm(n = num_states^2,
                                        # a = 0, b = 1,
                                        mean = 0,
                                        sd = 0.5),
         trade_positive = rtruncnorm(n = num_states^2,
                                     # a = 0, b = 1,
                                     mean = (1-ideal_point_distance), # negative of difference, implies trade more when ideal points closer
                                     sd = 0.5),
         trade_negative = rtruncnorm(n = num_states^2,
                                     # a = 0, b = 1,
                                     mean = ideal_point_distance, # positive of difference, implies trade less when ideal points closer
                                     sd = 0.5)) %>% 
  # Use logistic function to bound trade ties
  mutate_at(vars(trade_independent:trade_negative),
            ~logistic_function(x = ., k = 10, a = 0.75))


## Calculate dyadic trade flows as shares
climate_trade_shares <- climate_trade_data %>% 
  # Remove data for faux dyads -- where id_a == id_b
  mutate_at(vars(ideal_point_distance,
                 trade_independent, trade_positive, trade_negative),
            ~replace(., id_a == id_b, NA)) %>% 
  # Get each state's total trade
  group_by(id_a) %>% 
  mutate(sum_trade_independent = sum(trade_independent, na.rm = T),
         sum_trade_positive = sum(trade_positive, na.rm = T),
         sum_trade_negative = sum(trade_negative, na.rm = T)) %>% 
  # Standardize
  rowwise() %>% 
  mutate(share_dyadic_independent = trade_independent / sum_trade_independent,
         share_dyadic_positive = trade_positive / sum_trade_positive,
         share_dyadic_negative = trade_negative / sum_trade_negative) %>% 
  select(id_a, id_b, climate_ideal_a, climate_ideal_b, ideal_point_distance,
         share_dyadic_independent, share_dyadic_positive, share_dyadic_negative)
#' -- Interpretation
#' Each state has their own climate ideal point _a _b
#' The trade variables are:
#' How much _a exports to _b
#' share_dyadic_independent: what share of _a's exports are to _b


## Define a value of theta
theta_vector <- 0.75


## Calculate each state's utility to membership
climate_trade_utilities <- climate_trade_shares %>% 
  # Create a club's common climate goal
  mutate(theta = theta_vector,
         core_member_a = ifelse(climate_ideal_a >= theta, 1, 0),
         core_member_b = ifelse(climate_ideal_b >= theta, 1, 0)) %>% 
  # Flag trade to club members
  mutate(trade_to_club_independent = ifelse(core_member_b == 1, 
                                            share_dyadic_independent, NA),
         trade_to_club_positive = ifelse(core_member_b == 1, 
                                         share_dyadic_positive, NA),
         trade_to_club_negative = ifelse(core_member_b == 1, 
                                         share_dyadic_negative, NA),) %>% 
  # Get share of trade to all club members
  group_by(id_a) %>% 
  mutate(share_trade_club_independent = sum(trade_to_club_independent, na.rm = T),
         share_trade_club_positive = sum(trade_to_club_positive, na.rm = T),
         share_trade_club_negative = sum(trade_to_club_negative, na.rm = T)) %>% 
  ungroup() %>% 
  select(-starts_with("trade_to_club")) %>% 
  # Collapse the data to just _a's perspective, to calculate utilities
  select(id_a, climate_ideal_a, theta, core_member_a,
         starts_with("share_trade_club")) %>% 
  distinct() %>% # rows = 150
  # Calculate utilities to membership
  mutate(beta = 1) %>% # Define weights for climate versus trade
  rowwise() %>% 
  # Get most constrained core club member
  mutate(min_trade_exposure_independent = ifelse(
    core_member_a == 1, min(share_trade_club_independent), NA),
    min_trade_exposure_positive = ifelse(
      core_member_a == 1, min(share_trade_club_positive ), NA),
    min_trade_exposure_negative = ifelse(
      core_member_a == 1, min(share_trade_club_negative), NA)) %>% 
  ungroup() %>% 
  # Define that constraint as lambda
  mutate(lambda_independent = min(min_trade_exposure_independent, na.rm = T),
         lambda_positive = min(min_trade_exposure_positive, na.rm = T),
         lambda_negative = min(min_trade_exposure_negative, na.rm = T)) %>% 
  # Calculate utilities for the independent world
  # Utility to join
  mutate(utility_join_independent = (-abs(climate_ideal_a - theta)) + # Climate
           (1 - share_trade_club_independent) * (1 - lambda_independent)*beta + # Non-club trade
           (share_trade_club_independent*beta) - # Club trade
           (theta - lambda_independent)) %>%  # Aversion
  # Utility to abstain
  mutate(utility_abstain_independent = (-climate_ideal_a) + # Climate
           (1-share_trade_club_independent)*beta + # Non-club trade
           (share_trade_club_independent*(1-lambda_independent)*beta) + # Club trade
           (0)) %>%  # Aversion
  # Net, membership
  mutate(net_utility_independent = utility_join_independent - utility_abstain_independent,
         membership_independent = ifelse(net_utility_independent >= 0, 1, 0)) %>% 
  # Calculate utilities for the positive world
  # Utility to join
  mutate(utility_join_positive = (-abs(climate_ideal_a - theta)) + # Climate
           (1 - share_trade_club_positive) * (1 - lambda_positive)*beta + # Non-club trade
           (share_trade_club_positive*beta) - # Club trade
           (theta - lambda_positive)) %>%  # Aversion
  # Utility to abstain
  mutate(utility_abstain_positive = (-climate_ideal_a) + # Climate
           (1-share_trade_club_positive)*beta + # Non-club trade
           (share_trade_club_positive*(1-lambda_positive)*beta) + # Club trade
           (0)) %>%  # Aversion
  # Net, membership
  mutate(net_utility_positive = utility_join_positive - utility_abstain_positive,
         membership_positive = ifelse(net_utility_positive >= 0, 1, 0)) %>% 
  # Calculate utilities for the negative world
  # Utility to join
  mutate(utility_join_negative = (-abs(climate_ideal_a - theta)) + # Climate
           (1 - share_trade_club_negative) * (1 - lambda_negative)*beta + # Non-club trade
           (share_trade_club_negative*beta) - # Club trade
           (theta - lambda_negative)) %>%  # Aversion
  # Utility to abstain
  mutate(utility_abstain_negative = (-climate_ideal_a) + # Climate
           (1-share_trade_club_negative)*beta + # Non-club trade
           (share_trade_club_negative*(1-lambda_negative)*beta) + # Club trade
           (0)) %>%  # Aversion
  # Net, membership
  mutate(net_utility_negative = utility_join_negative - utility_abstain_negative,
         membership_negative = ifelse(net_utility_negative >= 0, 1, 0))

#### Analyze model outputs ####


## Count membership in each world
climate_trade_utilities %>% 
  select(starts_with("membership")) %>% 
  pivot_longer(names_to = "model", 
               values_to = "member", 
               membership_independent:membership_negative) %>% 
  group_by(model) %>% 
  summarize(members = sum(member)) %>% kable()

## Plot membership, independent
plot_independent <- climate_trade_utilities %>% 
  ggplot(., aes(climate_ideal_a, share_trade_club_independent, 
                shape = as.factor(membership_independent))) +
  geom_hline(data = climate_trade_utilities,
             aes(yintercept = lambda_independent), 
             color = "black", linetype = 2, linewidth = 2) +
  geom_vline(data = climate_trade_utilities,
             aes(xintercept = theta), 
             color = "#556B2F", linetype = 2, linewidth = 2) + 
  geom_point(color = "black", size = 5) +
  geom_point(aes(color = net_utility_independent), size = 4) +
  scale_colour_gradient2(low = "#5A1426", 
                         mid = "white",
                         high = "#556B2F") +
  scale_shape_manual(values = c(18, 16)) +
  labs(x = "Climate ideal point",
       y = "Export exposure to club, independent",
       color = "Net utility",
       shape = "Member") +
  guides(color=guide_colorbar(barwidth=20),
         shape=guide_legend(nrow=2,byrow=TRUE, order = 1)) +
  theme_clubs
plot_independent
ggsave(plot_independent, filename = "text/utilities_independent.pdf", units = 'in', height = 6, width = 8)


## Plot membership, positive
plot_positive <- climate_trade_utilities %>% 
  ggplot(., aes(climate_ideal_a, share_trade_club_positive, 
                shape = as.factor(membership_positive))) +
  geom_hline(data = climate_trade_utilities,
             aes(yintercept = lambda_positive), 
             color = "black", linetype = 2, linewidth = 2) +
  geom_vline(data = climate_trade_utilities,
             aes(xintercept = theta), 
             color = "#556B2F", linetype = 2, linewidth = 2) + 
  geom_point(color = "black", size = 5) +
  geom_point(aes(color = net_utility_positive), size = 4) +
  scale_colour_gradient2(low = "#5A1426", 
                         mid = "white",
                         high = "#556B2F") +
  scale_shape_manual(values = c(18, 16)) +
  labs(x = "Climate ideal point",
       y = "Export exposure to club, positive",
       color = "Net utility",
       shape = "Member") +
  guides(color=guide_colorbar(barwidth=20),
         shape=guide_legend(nrow=2,byrow=TRUE, order = 1)) +
  theme_clubs
plot_positive
ggsave(plot_positive, filename = "text/utilities_positive.pdf", units = 'in', height = 6, width = 8)

## Plot membership, negative
plot_negative <- climate_trade_utilities %>% 
  ggplot(., aes(climate_ideal_a, share_trade_club_negative, 
                shape = as.factor(membership_negative))) +
  geom_hline(data = climate_trade_utilities,
             aes(yintercept = lambda_negative), 
             color = "black", linetype = 2, linewidth = 2) +
  geom_vline(data = climate_trade_utilities,
             aes(xintercept = theta), 
             color = "#556B2F", linetype = 2, linewidth = 2) + 
  geom_point(color = "black", size = 5) +
  geom_point(aes(color = net_utility_negative), size = 4) +
  scale_colour_gradient2(low = "#5A1426", 
                         mid = "white",
                         high = "#556B2F") +
  scale_shape_manual(values = c(18, 16)) +
  labs(x = "Climate ideal point",
       y = "Export exposure to club, negative",
       color = "Net utility",
       shape = "Member") +
  guides(color=guide_colorbar(barwidth=20),
         shape=guide_legend(nrow=2,byrow=TRUE, order = 1)) +
  theme_clubs
plot_negative
ggsave(plot_negative, filename = "text/utilities_negative.pdf", units = 'in', height = 6, width = 8)


